Cite this Article
Zhang Jian, Tian Ming-Feng, Jin Guang-Xi, Xu Yuan-Feng, Dai Xi. Implementation of LDA+Gutzwiller with Newtons method. Chinese Physics B, 2017, 26(1): 017103  
Permissions
Implementation of LDA+Gutzwiller with Newtons method
Zhang Jian1, Tian Ming-Feng2, 3, Jin Guang-Xi4, Xu Yuan-Feng1, Dai Xi1, †
Beijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
LCP, Institute of Applied Physics and Computational Mathematics, Beijing 100088, China
CAEP Software Center for High Performance Numerical Simulation, Beijing 100088, China
Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621900, China

 

† Corresponding author. E-mail: daix@aphy.iphy.ac.cn

Abstract

In order to calculate the electronic structure of correlated materials, we propose implementation of the LDA+Gutzwiller method with Newton’s method. The self-consistence process, efficiency and convergence of calculation are improved dramatically by using Newton’s method with golden section search and other improvement approaches. We compare the calculated results by applying the previous linear mix method and Newton’s method. We have applied our code to study the electronic structure of several typical strong correlated materials, including SrVO3, LaCoO3, and La2O3Fe2Se2. Our results fit quite well with the previous studies.

PACS: 71.27+a;71.15.Mb;71.15.Nc;71.20.-b
Keyword:LDA+Gutzwiller;strongly correlated electrons;Newton’s method
1. Introduction

The density function theory (DFT) with local density approximation (LDA) have been successfully applied to calculate the band structures and related properties of many materials. However, the LDA fails to capture the strongly correlated nature of strongly correlated materials because it is based on a single particle picture. In order to solve this problem, several computational methods have been proposed, including , , and .

method has already been proposed and successfully utilized to investigate many important strongly correlated materials, such as ferromagnetic metals, iron superconductors, and topological Kondo insulators.[1,2,6,16] Compared to , is much faster and cheaper. However, the convergence of the self-consistent procedure is the bottle neck for this method, which prevents it from being applied to much bigger systems.

In this paper, we present an important improvement of the Gutzwiller method by combining it with Newton’s method (a numerical method to determine a function’s root) and golden section search (a technique for finding the extremum (minimum or maximum) of a unimodal function). Our hybrid method dramatically improves the convergence, and can make the calculations much faster and stable.

The rest of this paper is organized as follows. In Section 2, the theory behind the Gutzwiller method is introduced. In Section 3, implementation details are discussed, and the linear mix method and Newton’s method are compared. Section 4 details experimental results of our proposed method. Finally, Section 5 lists a short summary of our findings.

2. LDA+Gutzwiller derivation

In this section, we present the basic derivation of the Gutzwiller method. More detailed derivations and explanations can be found in our previous paper.[1,2]

2.1. Hubbard model

In this section, we derive the basic equations of the Gutzwiller approximation using the multi-orbital Hubbard model. The Hubbard model can be divided into two parts, one is the single particle ( ) part, which can be obtained from LDA calculations, the other is the interaction part ( ), which describes the electron–electron interaction by many-particle wave functions. Expanded expression is given as follows:

(1)
(2)
In Eq. (1), Greek letter subscripts denote a combined spin–orbital index, and the operator creates (annihilates) an electron on orbital α at site i. The first term of Eq. (1) describes electrons hopping between sites i and j, while the second term describes electrons hoping on one site between different orbitals. The last term describes the electron–electron interaction. In Eq. (2), m denotes strictly orbital index, while σ denotes spin index. The former three terms are density–density interaction terms: the first term describes interactions between two electrons in the same orbital but with different spins; the second term describes interactions between two electrons on different orbitals with different spins ( ); and the third term describes interactions between two electrons with the same spin but on different orbitals ( ). The third term’s factor is the smallest, thus electrons with the same spin tend to have lower system energy. This reflects Hund’s rules. The fourth term describes interactions of spin flipping and pair hopping.

In method, two types of interaction terms can be considered. One is called density–density Gutzwiller interaction (ddgw), in which only the first three terms in Eq. (2) are included. There are two parameters affecting the energy band, the quasi-particle weight R and the onsite energy revision Elm. R is between 0 and 1, and the band is compressed according to R. Elm affects the shift of the correlated bands. In this type, quasi-particle weight R and on-site energy Elm are vectors and Elm . α is the local orbital index. The other type, rotational invariable Gutzwiller interaction (rtgw), includes all interacting terms. In rtgw, quasi-particle weight and on-site energy are matrix and Elm . α and β are local orbital index. In the present paper, we will focus on the density–density type just for simplicity. The generalization to the rotational invariant type is straight forward. The result of the rotational invariant type is shown in Section 4.2.

2.2. Gutzwiller variational wave function

The Gutzwiller wave function (GWF) [3] is a many-particle wave function. We act the Gutzwiller projection operator upon the Hartree uncorrelated wave function (HWF)

(3)
This operator P is composed by projector P i
(4)
where i is the site index, and P i projects the wave function to many-particle configurations
(5)

Here, I denotes the Fock states, the parameters c I are variational parameters, which adjust the weight of each configuration in the wave function. The GWF is obtained by setting the proper parameters c I .

2.3. Variational problem constraints

There are two variational problem constraints. First, the GWF must be normalized as follows:

(6)
The other constraint is[3]
(7)
and we denote
(8)
(9)

2.4. Lagrange equation

The Lagrange function is , in which c I is the Lagrange variation parameter to determine the GWF. is the Lagrange variation parameter to determine effective one-particle wave function . is the eigen wave function from LDA calculation, we use as single particle basis. and are the Lagrange multipliers of the normalization constraint for the effective one-particle wave function and the GWF, respectively. is the Lagrange multiplier of the particle number operator constraint, and its physical meaning is each orbit’s chemical potential.

(10)

In this non-linear equation, both the effective single particle wave function and the many-particle wave function c I are unknown, and they are related to each other. To solve this equation, we designed a two-step self-consistent loop. First, all parameters related to the many-particle wave function ( , Elm ) are set as constants, and we can find a proper effective single particle wave function . Second, all parameters related to the single particle wave function ( , ) are set as constants, and we can constitute hamiltonian on the non-interacting many-particle wave function and find the proper variational parameters c I . The equation for the effective single-particle part is as follows:

(11)

The equation for the many-particle part is

(12)

3. Implementation and optimization

In this section, we focus on implementation details and optimization techniques applied to the Gutzwiller method. We introduce the interface of method. Then outline the two self-consistent loops (outer loop and inner loop), and compare the linear mixed method and Newton’s method. Lastly, we discuss the golden section search method used in the inner loop.

3.1. Interface of LDA+Gutzwiller

The Gutzwiller method requires two input files. One is the eigenvalue of each k point E nk , which results from LDA calculation. The other is the overlap matrix , . is the eigenfunction from LDA calculation. is the local orbital. In our method, strong correlation confined to each site is required, therefore local orbitals between neighbor sites should not be overlapped. The local orbital can be a Wannier function[4,5] of the LDA bloch wave, or atomic local orbitals.

Tight binding Harmiltonian can be built on local orbital basis. Since some tight binding models are provided in real space, the Hamiltonian is , and it is better to transform to reciprocal space by Fourier transformation, to facilitate the effective single particle calculation. Details on building local orbital and projection can be found in previous paper.[6]

3.2. Two-step loop of Gutzwiller method

In the sigle-particle part (Section 3.3), with given input and Elm , we can find proper . Then we can calculate parameters related with , such as occupation number and derivative . In the many-particle part (Section 3.4), with occupation number and derivative from single-particle part, we can find proper c I . The inner self-consistent loop is only related to many-particle part (Eq. (12)). In the inner loop, with given occupation number and derivative , we need to find proper chemical potential of each orbit to meet the constraint condition (Eq. (7)). Then we use occupation number , derivative and chemical potential (from inner self-consistent loop) to find proper c I . Then we can calculate the output and Elm . The two-step loop is shown in Fig. 1.

Fig. 1. (color online) Flow chart of the two-step self consistent loop.
3.3. Effective single-particle equation

From Eq. (11), we get an effective single particle Hamiltonian

(13)
(14)
(15)
is quasi-particle weight, Elm is the chemical potential of local orbital α. Both are functions of c I . is the occupation number of local orbital ( ), is the kinetic energy, and f nk is the Fermi dispersion of the effective single particle wave function. is either one or zero, which denotes eliminating an electron on orbital α in Fock state , and then mapping to Fock state I .

3.4. Many-particle equation

From Eq. (12), we obtain a many-particle Hamiltonian written on the basis of Fock states,

(16)
(17)
(18)
where
(19)
and are two input parameters gotten from the effective single particle equation. denotes the electronic occupation on orbital α in the Fock state I ( ). The ground state eigen vector of H m is the Lagrange variational parameter c I . Then we can calculate ( ) to see whether the constraint has been met ( ). Then we adjust to meet the constraint condition given by Eq. (7).

3.5. Outer self-consistent loop and Newton’s method

Each self-consistent loop was separated into two parts. In the effective single-particle equation, input is quasi-particle weight and on-site energy Elm , output is kinetic energy and effective occupation number . Conversely, in the many-particle equation, input is kinetic energy and occupation number , output is quasi-particle weight and on-site energy Elm . The two parts constitute a complete self-consistent loop. We can view the loop as function of , with input and , and function results of and . For convenience, input parameters are called , while the output parameters are called

(20)
This equation is simplified as
(21)
Previously, we use linear mix method, that is,
(22)
where λ is a real number between 0 and 1.

To make this self-consistent procedure more efficient, we then applied Newton’s method. Using this method, the self-consistent problem (function g) is changed to a root finding problem (function f) and defined as

(23)
The self-consistent point is the point . Taking the derivative at point , which is near the root point , we get
(24)

Taking , becomes the next input point for our consistent loop. This equation becomes

(25)
and the next input point is
(26)
is the Jacobi matrix. To apply Newton’s method, we must derive the Jacobian: , , and . In order to get this Jacobi matrix, we break it into two parts: the effective single particle part , , , (see Appendix A), and the many-particle part , , , (see Appendix B).

In the many-particle part, there is a intermediate variable , and we need and to yield the Jacobi matrix for the many-particle part. The details for these derivatives are shown in Appendix C.

3.6. Inner self-consistent loop and Newton’s method

The inner self-consistent loop corresponds to the many-particle part outlined in Section 3.4. The purpose is to find the proper in Eq. (16). In this loop, the parameter of many-particle Hamiltonian H m is the chemical potential , , and are given by single-particle part (Eq. (11)). With given , we can get Hamiltonian H m . Diagonalizing this Hamiltonian yields the ground state wave function C I and the corresponding , . We define a function , where input is and output is , .

To meet the constraint , we define a new function

(27)

The Jacobi matrix for the inner loop is . For derivation details, please refer to Appendix C and Eq. (A10).

3.7. Golden section search with Newton’s method

We also apply golden section search[24] to improve the convergence of our code. Newton’s method sets direction in the parameter space, and then golden section search sets the step length. Newton’s method is very sensitive to the starting point. In some situations, it can cost more steps or cannot find solution, if Newton’s method starts with a bad point. With the technique of golden section search, in each step, we can find the point nearest to the solution on the step’s Newton’s line, which is a better starting point for next step, thus it can improve the convergence. Results are shown in Section 4.3.

4. Example

In this section, we show results from SrVO3, LaCoO3 and La2O3Fe2Se2, and compare the LDA and LDA+G band structures. We will also compare our method to the traditional Gutzwiller method. Convergence behavior will be shown with respect to convergence steps.

4.1. SrVO3

SrVO3 is a 3d orbital correlated metal. Its simple cubic perovskite structure and non-magnetic electronic state make it an ideal test material.[5,7,8] The LDA result is inconsistent with experiments.[9–11] These features can be understood through our calculations. We choose a Wannier function as the interacting local orbits with interaction energy U = 5 eV, . Figure 2 shows the band structure calculated by . The band width using our method is reduced by . These results are in good agreement with previous experiments.[9–11]

Fig. 2. (color online) The band structure of SrVO3 from the LDA and the result. U = 5 eV and .

The results of quasi-particle weight and on-site energy by our new application are all well in comparable with the traditional application. Figure 3 shows convergence behavior of the two methods in case of SrVO3. Case far away from the transition point ( , ) and near the transition point ( , eV) are shown in upper and lower part respectively. In both cases, applying Newton’s method reduces the number of convergence steps.

Fig. 3. (color online) The Gutzwiller energy of SrVO3 with respect to steps. Upper case: U = 1, J = 0. Lower case: U = 6, J = 0. Figures on the left are zoomed images of right figures. Figures on the left are zoomed images of the right figures.

SrVO3 showed one electron in the 3d orbitals, consistent with previous findings.[13] The Brinkman–Rice transition[14] occurred at certain combination of U and J. As shown in Fig. 4, the Gutzwiller method can converge into two different results, which is also consistent with previous research.[15] The previous linear mix application cannot find self-consistent solution for branches with small quasi-particle weight. The application of Newton’s method can improve the convergence of the Gutzwiller method.

Fig. 4. The quasi-particle weight R[12] with respect to U. From right to left, , 0.1, 0.0.
4.2. LaCoO3

LaCoO3 is a typical correlation material. Since there are two strongly correlated atoms in one cell, and each atom has 10 local orbitals, convergence of this case is more difficult than SrVO3. We calculated bulk LaCoO3[16] band structure with LDA and LDA+G (rtgw). We applied rotational invariant code for this case. The Gutzwiller method adjusted crystal field splitting between and due to the effect of onsite energy Elm, and renormalized and bandwidth using quasi-particle weight R. LDA and LDA+G bands are plotted in Fig. 5. The Gutzwiller method shows that LaCoO3 is semiconductor, while LDA calculation suggests it is metal. The band gap was approximately 0.36 eV, which is in good agreement with experiments.[17–19]

Fig. 5. (color online) LaCoO3 band structure with interaction , Hund’s rule coupling .
4.3. La2O3Fe2Se2

La2O3Fe2Se2 is considered a strongly correlated material and is claimed to be a semiconductor and Mott insulator.[20] We consider the magnetic order of La2O3Fe2Se2. Each cell has 16 strongly correlate atoms, and each atom has 10 local orbitals. This makes convergence very difficult. The band structure with interaction and is shown in Fig. 6. The band gap was approximately 0.121 eV, which agrees with past findings.[21–23]

Fig. 6. (color online) La2O3Fe2Se2 band structure with interaction , Hund’s rule coupling .

Calculations were performed using both linear mix and Newton’s methods. The linear mix method cannot achieve convergence. In the inner loop, the linear mix method cannot locate appropriate chemical potentials for each orbital defined by Eq. (12) with the constraint condition given by Eq. (7). In contrast, Newton’s method locates appropriate chemical potentials within 6 steps. Figure 7 shows the average absolute deviation of the difference between and with respect to number of steps. Figure 8 shows the Gutzwiller energy in the inner loop with respect to the number of steps.

Fig. 7. (color online) The average absolute deviation of the difference between and with respect to number of steps.
Fig. 8. (color online) The Gutzwiller total energies in the inner self-consistent loop with respect to inner loop steps, implementation of linear mix and Newton’s method respectively.
5. Conclusion

Our implementation of LDA+G with Newton’s method, golden section search, and other revision of previous Gutzwiller implementation were detailedly described in this paper. First, the Hubbard model (Section 2.1) and Gutzwiller wave function (Section 2.2) were introduced. Then, we put forward a two-step self-consistent method to solve the nonlinear Lagrange equation (Section 2.4), and apply Newton’s method to both the outer and inner loops (Sections 3.5 and 3.6 respectively). Finally, we tested our method on SrVO3, LaCoO3, and La2O3Fe2Se2. Results for all materials agree with both the experimental data and previous calculation. Conclusively, this proposed method dramatically improves the efficiency and convergence of the implementation of Gutzwiller method.

Reference
[1] Deng X Y Dai X Fang Z 2008 EPL 83 37008
[2] Deng X Y Wang L Dai X Fang Z 2009 Chin. Phys. B 79 075114
[3] Bünemann J Weber W Gebhard F 1998 Phys. Rev. B 57 6896
[4] Anisimov V I Kondakov D E Kozhevnikov A V Nekrasov I A Pchelkina Z V Allen J W Mo S K Kim H D Metcalf P Suga S Sekiyama A Keller G Leonov I Ren X Vollhardt D 2005 Phys. Rev. B 71 125119
[5] Lechermann F Georges A Poteryaev A Biermann S Posternak M Yamasaki A Andersen O K 1998 Phys. Rev. B 74 125120
[6] Zhao J Z Zhuang J N Deng X Y Yan Bi Cai L C Fang Z Dai X 2012 Chin. Phys. B 21 057106
[7] Liebsch A 2003 Phys. Rev. Lett. 90 096401
[8] Nekrasov I A Keller G Kondakov D E Kozhevnikov A V Th Pruschke Held K Vollhardt D Anisimov V I 2005 Phys. Rev. B 72 155106
[9] Fujimori A Hase I Namatame H Fujishima Y Tokura Y Eisaki H Uchida S Takegahara K de Groot F M F 1992 Phys. Rev. Lett. 69 1796
[10] Inoue I H Hase I Aiura Y Fujimori A Haruyama Y Maruyama T Nishihara Y 1995 Phys. Rev. Lett. 74 2539
[11] Sekiyama A Fujiwara H Imada S Suga S Eisaki H Uchida S I Takegahara K Harima H Saitoh Y Nekrasov I A 2004 Phys. Rev. Lett. 93 156402
[12] Bünemann J Gebhard F Thul R 2003 Phys. Rev. B 67 075103
[13] Medici Luca de’ 2011 Phys. Rev. B 83 205112
[14] Brinkman W F Rice T M 1907 Phys. Rev. B 2 4302
[15] Lanatà N Strand H U R Dai X Hellsing B 2012 Phys. Rev. B 85 035133
[16] Wang Y L Wang Z J Fang Z Dai X 2015 Phys. Rev. B 91 125139
[17] Chainani A Mathew M Sarma D D 1992 Phys. Rev. B 46 9976
[18] Abbate M Fuggle J C Fujimori A Tjeng L H Chen C T Potze R Sawatzky G A Eisaki H Uchida S 1993 Phys. Rev. B 47 16124
[19] Arima T Tokura Y Torrance J B 1993 Phys. Rev. B 48 17006
[20] Free D G Evans J S O 2010 Phys. Rev. B 81 214433
[21] Zhu J X Yu R Wang H d Zhao L L Jones M D Dai J H Abrahams E Morosan E Fang M H Si Q M 2010 Phys. Rev. Lett. 104 216405
[22] Lei H C Bozin E S Llobet A Ivanovski V Koteski V Belosevic C J Cekic B Petrovic C 2012 Phys. Rev. B 86 125122
[23] Jin G X Wang Y L Dai X He L X 2016 Phys. Rev. B 94 075150
[24] Press W H Teukolsky S A Vetterling W T Flannery B P 2007 Numerical Recipes 3 The Art of Scientific Computing New York Cambridge University Press Section 10.2